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Abstract 


An  explicit,  second-order  accurate,  total  variation 
diminishing  (TVD)  scheme  and  the  MacCormack  scheme  are  applied 
to  the  Euler  equations  in  axisymmetric  form  to  study 
hypersonic  blunt-body  flows.  The  modif ied-f lux  approach  of 
Harten,  with  modification  by  Yee1,  for  two-dimensional  flows 
is  extended  to  treat  axisymmetric  flows.  Calculated  flow 
properties  for  the  steady-state,  blunt-body  problem  such  as 
shock  standoff  distance,  bow-shock  shape,  surface  pressure 
distribution  and  entropy  jump  conditions  are  compared  to 
theory,  results  of  the  MacCormack  scheme,  and  experimental 
data  for  Mach  numbers  of  3.C,  4.03,  5.06,  6.03  and  8.1. 
Additionally,  the  TVD  and  MacCormack  schemes  are  used  to 
simulate  numerically  the  unsteady  shock  impingement  on  a 
sphere.  Results  are  compared  to  experimental  data  for  a  shock 
Mach  number  of  2.89.  Analysis  of  the  numerical  simulations 
provide  suitable  ranges  of  values  for  the  entropy  correction 
parameter  and  the  Courant  ( CFL)  number.  A  brief  comparison  of 
limiters  for  the  unsteady  problem  is  also  presented.  The 
high-resolution,  shock-capturing  capability  and  robustness  of 
the  TVD  scheme  is  clearly  shown. 


ANALYSIS  OF  HYPERSONIC  BLUNT-BODY  FLOWS 

USING  A 

TOTAL  VARIATION  DIMINISHING  (TVD)  SCHEME 

AND  THE 

MACCORMACK  SCHEME 


I •  Introduction 

Programs  of  national  significance  such  as  the  National 
Aerospace  Plane  (NASP) ,  the  Strategic  Defense  Initiative  (SDI) 
and  the  Aeroassisted  Orbital  Transfer  Vehicle  (AOTV)  have 
renewed  interest  in  hypersonic  aerodynamics.  Hypersonic  wind 
tunnels  and  ground-test  facilities  developed  in  the  1960 's 
during  the  initial  era  of  hypersonics  research  interest  fell 
into  disuse  following  the  drastic  funding  cuts  of  the  1970's 
and  are  just  now  being  rebuilt  and  refurbished2.  Due  to  these 
constraints  on  physical  testing  of  hypersonic  vehicle  models, 
numerical  simulations  become  particularly  important. 
Computational  fluid  dynamics  (CFD)  offers  immediate  means  for 
hypersonic  aerodynamics  research,  development  and  design  for 
a  number  of  reasons.  There  are  three  main  advantages  of  CFD 
over  field  tests  and  laboratory  experiments3: 

1)  lower  cost, 

2)  increased  data  collection  capability  -  virtually 
any  flow  quantity  can  be  computed  in  the  field  of 
interest,  and 

3)  exact  repeatability,  precise  control  over  initial 
and  boundary  conditions,  and  the  related  capability 
to  independently  vary  important  parameters. 
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Blunt-body  flows  are  of  particular  interest  because 
hypersonic  vehicles  are  designed  with  blunt  noses  to  reduce 
aerodynamic  heating.  The  solution  of  the  Euler  equations  of 
fluid  flow  in  the  nose  region  of  a  blunt-body,  such  as  those 
presented  in  this  study,  provide  numerical  estimates  of  shock 
standoff  distance,  bow-shock  shape,  and  surface  pressure 
distributions.  Additionally,  these  solutions  provide  edge 
properties  for  input  to  boundary  layer  codes,  which  in  turn 
can  provide  estimates  of  heat  transfer  and  skin  friction.  In 
the  study  of  transient  flow  phenomena,  application  of  the 
numerical  schemes  to  the  unsteady  shock  impingement  problem 
can  provide  pressure  and  density  contours  about  a  sphere  as  a 
blast-wave  passes  over  it.  This  ability  to  accurately 
simulate  the  complex  shock-diffraction  process  is  particularly 
important  for  the  analysis  of  a  building's  ability  to  maintain 
structural  integrity  during  the  passage  of  a  blast-wave  from 
a  nuclear  explosion. 

1 . 1  TVD  and  MacCormack  Scheme  Background 

In  this  section,  the  background  information  is  condensed 
from  the  excellent  survey  reported  by  Yee  in  Reference  1. 

For  complex  flowfields,  monotone  and  first-order  upwind 
schemes  are  too  diffusive.  Monotone  schemes  produce  smooth 
transitions  near  discontinuities,  but  they  are  only  first- 
order  accurate.  High-resolution  shock-capturing  numerical 
simulations  are  not  possible  with  these  schemes  on  grids  of 
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reasonable  size.  Schemes  of  higher-order  accuracy  have  been 
developed  over  the  past  decade  and  the  development  of  such  is 
an  area  of  active  research. 

Classical  shock-capturing  schemes  use  linear  numerical 
dissipation  -  the  same  amount  everywhere  or  consists  of 
adjustable  parameters  for  each  problem.  Modern  schemes  employ 
numerical  dissipation  in  a  nonlinear  fashion  -  the  amount 
varies  from  grid  point  to  grid  point  and  is  implemented 
automatically  within  the  scheme  with  few,  if  any,  adjustable 
parameters1.  This  technique  was  developed  to  overcome  the 
following  deficiencies: 

1)  development  of  spurious  oscillations  whenever  the 
solution  is  not  smooth, 

2)  development  of  nonlinear  instabilities  when 
discontinuities  are  encountered,  and 

3)  the  selection  of  a  nonphysical  solution. 

Total  variation  diminishing  (TVD)  schemes  are  a  class  of 
modern,  high-resolution,  shock-capturing  schemes  introduced  by 
Harten4.  A  main  objective  of  this  class  of  schemes  is  to 
ensure  that  solutions  of  the  Euler  equations  converge  to 
physically  correct  solutions.  This  is  extremely  important  if 
one  is  dealing  with  fluid  flows  for  which  analytic  solutions 
or  experimental  data  are  not  readily  available.  The  basic 
idea  behind  the  notion  of  a  TVD  property  is  that  the  total 
variation  of  the  numerical  solution  of  an  initial-value 
problem  will  not  grow  as  the  solution  evolves  in  a  time-like 
fashion.  See  Reference  4  for  precise  mathematical  details. 
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It  should  be  emphasized  that  the  TVD  property  is  only  valid 
for  systems  of  homogeneous  scalar  hyperbolic  equations.  For 
nonhomogeneous  hyperbolic  equations,  in  order  for  the  source 
term  to  not  influence  the  TVD  property,  it  is  restricted  to  a 
special  class  of  functions  and  fluid  flows1.  For  finite- 
volume,  upwind  TVD  schemes,  similar  to  the  one  employed  in 
this  study,  Wang  and  Widhopf5  prove  the  TVD  property  is 
satisfied  for  the  axisymmetric  version  of  the  Euler  equations 
written  in  conservation  form. 

The  MacCormack  scheme  has  been  an  industry  standard  for 
many  years  and  details  of  the  numerical  algorithm  other  than 
boundary  conditions  are  not  reported  in  this  study.  The 
scheme  has  been  applied  to  a  wide  variety  of  fluid  flow 
problems.  Application  of  the  scheme  to  the  Navier-Stokes 
equations  for  axisymmetric  flows  with  finite-rate  chemical 
kinetics  has  been  reported  by  Shang  and  Josyula6.  The 
MacCormack  scheme  applied  to  the  Euler  equations  in 
axisymmetric  form  was  used  for  this  study. 

As  will  be  shown  in  the  sections  forthcoming, 
comparisons  between  flowfield  properties  computed  with  the  TVD 
scheme  and  the  MacCormack  scheme  are  very  good.  Comparisons 
between  the  schemes,  theory,  and  experiment  are  also  very 
good. 
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1 . 2  Problem  Statement 


The  main  objectives  of  this  study  are  to  show  the  high- 
resolution  and  robustness  capabilicies  of  the  two-dimensional 
TVD  scheme  extended  to  axisymmetric  form.  Applications  to  the 
steady-state,  blunt-body  problem  and  the  problem  of  unsteady 
shock  impingement  on  a  sphere  provide  a  wide  range  of 
flowfield  features  possible  with  the  Euler  equations. 
Furthermore,  within  the  context  of  each  application,  suitable 
values  of  parameters  inherent  to  using  the  scheme,  such  as 
Courant  (CFL)  number,  entropy  correction  parameter,  8p  , 
appropriate  limiters  as  well  as  boundary  and  initial 
conditions,  are  numerically  investigated. 

For  the  steady-state,  blunt-body  problem,  flowfield 
quantities  computed  with  the  TVD  scheme  will  be  compared  to 
theory,  experimental  data,  and  results  of  the  MacCormack 
scheme.  Specific  items  considered  are  the  shock  standoff 
distance,  bow-shock  shape,  body  surface  pressure  distribution, 
and  the  entropy  jump  condition. 

For  the  unsteady  impingement  of  a  shock  on  a  sphere, 
flowfield  quantities  computed  with  the  TVD  scheme  will  be 
compared  to  theory  and  experimental  data.  Specific  items 
considered  are  the  density  and  pressure  distributions  in  the 
flowfield  and  the  resolution  of  the  complex  interactions. 
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II .  Analysis 


A  discussion  of  the  governing  equations,  the 
transformation  from  the  physical  domain  to  the  computational 
domain,  boundary  conditions,  initial  conditions,  theoretical 
predictions,  and  a  description  of  the  convergence  criteria 
will  be  presented  in  the  next  few  sections. 

2 . 1  Governing  Equations 

The  Euler  equations  are  statements  of  conservation  of 
mass,  momentum  and  energy.  These  principles  are7: 

1)  mass  can  be  neither  created  nor  destroyed, 

2)  the  time  rate  of  change  of  momentum  of  a  body 
equals  the  net  force  exerted  on  it,  and 

3)  energy  can  be  neither  created  nor  destroyed;  it  can 
only  change  in  form. 

The  properties  of  the  fluid  medium  are  characterized  by 
density  p,  the  pressure  p,  the  internal  energy  per  unit 
mass  e,  and  the  fluid  velocity  u  as  functions  of  position  and 
time.  An  equation  of  state  is  assumed  of  the 

form  e  =  f(p,  p)  .  Additionally,  it  is  assumed  that  no  viscous 
forces,  bod>  forces,  heat  conduction  or  energy  sources  are 
present.  The  conservation  form  of  the  two-dimensional  Euler 
equations  can  be  given  in  conservative  variables 
p,  pu,  pv,  and  Et.  The  equations  have  the  ability  to  describe 
internal  discontinuities  such  as  shock  waves 
(where  p,  pu,  pv,  and  e  are  discontinuous)  and  contact 
surfaces  (where  p,  and  e  are  discontinuous)  which  are 
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frequently  encountered  in  fluid  flow8,9.  The  Euler  equations 
can  be  written  in  flux  vector  form: 


dU  +  dE(U)  +  dF{U) 
dt  dx  dy 


where 


pyP' 

mxy* 

mxyP 

(p  +  ml/ p )  yp 

™yy* 

E  = 

/nxvyP 

Ety\ 

mjp  ( Et  +  p)  yP 

myy* 

myuy* 

(p  +  my/ p )  yP 
wy/p  ( Et  +  p)yp 


0 

0 


0 


(2) 


and  P  =  0,1  for  two-dimensional  and  axisymmetric  cases, 

respectively;  mx  =  p u,  my  =  pv,  the  total  energy,  Et,  is 


E 


t 


pe  + 


ml  +  ml 

2p 


(3) 


and  pressure  is  given  as 

p  =  (y-l)  pe 


(4) 
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2 . 2  Coordinate  Transformation 


For  the  numerical  solution  of  these  partial  differential 
equations  (PDEs) ,  it  is  imperative  that  boundary  conditions  be 
accurately  represented.  In  fact,  the  boundary  conditions  are 
what  differentiate  one  fluid  flow  problem  from  another,  as  far 
as  the  governing  equations  are  concerned10.  In  general,  the 
Cartesian  coordinate  system  is  limited  in  usefulness  to 
configurations  that  can  be  represented  by  a  rectilinear  shape. 
For  a  more  general  body  shape,  a  transformation  of  the 
governing  equations  to  a  curvilinear,  boundary-conforming 
system  is  more  appropriate.  A  transformation  of  the 
form  i  =  £(x,y)  and  tj  =  r\{x,y)  is  used  to  map  Eq  (1)  from  a 
general  physical  domain  ( x,y )  to  a  rectangular  computational 
domain  (£,ti)  as  illustrated  in  Figure  1.  The  transformed 
strong  conservation  form  of  the  Euler  equations,  assuming  the 
computational  grid  is  fixed  in  time,  is 

dU  +  dE(U)  +  dF(U)  =  g  (5) 

d t  dx  dy 

where 


U  =  JU 

E  =  JiixE  +  5^) 
F  =  +  r\yF) 

S  =  Ss 


and  J  =  x^yA  -  x^y^  is  the  Jacobian  of  the  transformation.  The 
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Figure  1.  Mapping  from  Physical  to  Computational  Domain 


solution  via  Roe's  approximate  Riemann  solver11  (sometimes 
referred  to  as  the  "local  characteristic  approach")  requires 
evaluation  of  the  Jacobians  A  and  B  of  E  and  F.  Through  the 
application  of  the  chain  rule  these  can  be  written  as 

*-<5*A  +  5yB>  (7) 

B  =  (t)^  +  n  yB) 


where  A  =  dE/dU  and  B  =  dF/dU.  The  eigenvalues  of  A  and  B  are 


a 


h 


!hxu  +  hyv  -  khc' 
hxu  +  hyv 
hxu  +  hyv  +  khc 


hxu  +  hyv 


) 


where 


h  =  £  for  / 
h  =  t)  for  B 


and 


hxu  +  =  Uc 

C  =  Hp/  p 

**  =  ^  +  hY 


(8) 


(9) 


(10) 


where  Uc  is  referred  to  as  the  contravariant  velocity.  The 
contravariant  base  vectors  of  the  curvilinear  coordinate 
system  are  normal  to  the  coordinate  surfaces12,  and  the 
contravariant  velocities  are  the  inner  product  of  the  base 
vectors  with  the  velocity  vector. 
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The  right  eigenvectors  of  A  and  B,  (R^,  Rl,  r\,  R^)  ,  are 


1.0 

1.0 

u  -  khlc 

u 

“  : 

V  -  kh2C 

II 

<N  XJ 
* 

V 

H  -  khluc  -  kh2vc 

—  ( u2  +  v 2) 

2  j 

1.0 

0.0 

q 

u  +  lcA1c 

A 

~kh2 

Rh  = 

V  +  kh2C 

Rh  = 

khi 

H  +  fcA1uc  + 

khlV~kh2U_ 

where 


H 

khi 

kh2 


VP  + 
(Y-l) p 

K 

Jh2x  +  h2y 

K 

JR  *  R 


+  v2) 


(12) 


Furthermore,  the  total  enthalpy  and  the  speed  of  sound  can  be 
written  in  terms  of  the  conserved  variables  as 

H  -  (Et  +  — ) 

P 

p  =  (y-l)  (Et  -  {ml  +  w2)  /  (2p)  )  (13) 

c2  =  (y-l)  (//  -  (/??*  +  /i?y)  /  (2p2)  ) 

For  the  actual  implementation  of  the  equations,  the  flow 
quantities  were  nondimensionalized  in  a  consistent  fashion  as 
shown  in  the  Appendix. 
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2 . 3  Boundary  Conditions 


The  finite-volume  TVD  scheme  constructs  cells  using  node 
points  as  vertices  in  the  physical  domain.  The  ensuing  grid 
of  cell  centers  is  then  used  for  actual  calculations. 
Furthermore,  mirror  images  (referred  to  as  ghost  points  in 
this  study,  sometimes  also  referred  to  as  phantom  points)  are 
placed  at  cell  centers  just  outside  the  physical  domain  to  aid 
in  the  enforcement  of  the  boundary  conditions.  The  finite- 
difference  MacCormack  scheme,  however,  uses  actual  node  points 
for  calculations  and  does  not  use  ghost  points. 

A  general  mapping  of  the  physical  domain  to  the 
computational  domain  was  illustrated  in  Figure  1.  The 
physical  domain  shown  was  that  for  the  steady-state 
calculations  whereby  a  blunt-body  with  a  spherical  nose  is 
modelled.  Referring  to  Figure  1,  the  physical  domain  consists 
of : 

1)  a  line  of  axisymmetry  (front  stagnation  line), 

2)  a  parabolic  outer  boundary  where  freestream  conditions 
are  imposed, 

3)  a  supersonic/hypersonic  outflow  boundary,  and 

4)  the  spherical  body  surface. 

The  parabolic  shape  of  the  outer  boundary  was  selected  based 
on  the  fact  that  bow-shock  shapes  are  parabolic.  The  intent 
was  to  minimize  the  number  of  grid  points  between  Lhe 
freestream  boundary  and  the  captured  shock. 
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For  the  unsteady  calculations,  the  physical  domain  is  a 
half-circle  (sphere)  as  will  be  shown  in  Section  VI.  The 
physical  domain  consists  of: 

1)  front  and 

2)  rear  lines  of  axisymmetry  (stagnation  lines) , 

3)  a  circular  outer  boundary  where  freestream  conditions 
are  imposed,  and 

4)  the  spherical  body  surface. 

The  boundary  conditions  are  described  in  detail  in  the  next 
few  sections. 

2.3.1.  Line  of  Axisvmmetrv.  For  the  TVD  scheme,  the 
reflection  principle  was  applied  across  the  line  of 
axisymmetry.  Derivatives  of  conserved  variables  normal  to  the 
line  of  axisymmetry  and  the  v  com^nent  of  velocity  were 
specified  to  vanish.  Furthermore,  since  the  physical  area  on 


the  line  of 

axisymmetry 

is 

zero,  the 

flux 

at 

that  cell 

interface  was  also  set 

to 

zero.  Exactly 

how 

this 

is 

implemented 

in  the  numerical 

algorithm 

will 

be 

shown 

in 

Section  IV.  For  the  MacCormack  scheme,  a  limiting  form  of  the 
governing  axi symmetric  equations,  developed  by  Shang  and 
Josyula6,  was  used  to  overcome  the  geometric  singularity  at 
the  line  of  axisymmetry  encountered  in  the  finite-difference 
formulation.  As  will  be  shown  in  the  steady-state  solutions 
of  the  Euler  equations  in  Section  V,  oscillations  of  flowfield 
quantities  in  the  stagnation  region  near  the  line  of 
axisymmetry  are  a  common  occurrence.  In  some  cases,  these 
oscillations  lead  to  numerical  instabilities. 
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2.3.2.  Freestream.  It  is  assumed  that  the  parabolic 
outer  boundary  for  the  steady-state,  blunt-body  problem,  and 
the  circular  outer  boundary  for  the  unsteady  shock  impingement 
problem,  are  located  sufficiently  far  away  from  the  body 
surface  that  freestream  boundary  conditions  can  be  Imposed. 
For  the  unsteady  calculations,  it  was  necessary  to  track  the 
planar  moving  shock  wave  and  impose  pre-  and  post-shock 
boundary  conditions  consistent  with  the  Rankine-Hugoniot 
relations  for  changes  across  a  normal  shock. 

2.3.3.  Outflow.  The  outflow  boundary  condition  described 
in  this  section  pertains  to  the  steady-state,  blunt-body 
problem  shown  in  Figure  1.  At  the  top  of  the  sphere,  the 
fluid  flow  has  accelerated  to  a  supersonic  Mach  number  from  a 
subsonic  Mach  number  in  the  nose  region  immediately  behind  the 
normal  shock.  Two-point  extrapolation  was  used  for  the 
outflow  boundary  values  in  the  TVD  scheme.  A  "no-change"  or 
one-point  extrapolation  method  worked  just  as  well  and  was 
used  in  the  MacCormack  scheme.  One  advantage  of  ^wo-point 
extrapolation  is  that  contours  of  flowfield  values  are 
smoother  at  the  outflow  boundary  for  relatively  coarse  grids. 
This  will  be  apparent  in  the  Mach  contour  plots  to  be  shown  in 
Section  V. 
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2.3.4.  Body  Surface.  The  following  paragraph  is  quoted 


directly  from  Anderson  et  al8: 

We  now  address  the  problem  of  applying  surface  boundary 
conditions  when  shock-capturing  methods  are  used.  As 
Moretti  (1969)  has  pointed  out,  the  correct  application 
of  boundary  conditions  is  a  difficult  task.  Incorrect 
wall  conditions  can  provide  locally  polluted  results 
and  in  many  cases  can  destroy  a  solution.  Hyperbolic 
equations  are  particularly  sensitive.  Due  to  their 
wave-like  nature,  boundary  errors  are  propagated 
throughout  the  mesh  reflecting  until  actual  instability 
can  result. 

For  these  reasons,  particular  care  was  taken  when  selecting 
the  following  body  surface  boundary  conditions. 

At  the  sphere  surface,  a  "no-flow  through"  condition  was 
imposed.  That  is,  the  normal  component  of  velocity  is  set  to 
zero  while  the  tangential  component  is  preserved.  For  the  TVD 
scheme,  ghost  points  inside  the  physical  body  surface  were 
used  and  these  conditions  are 


u 

u 


nl 

cl 


-U 


n2 


U 


t2 


(14) 


where  j  =  1  corresponds  to  the  r|  coordinate  in  the  radial 
direction  for  the  ghost  point  as  shown  in  Figure  2.  The 
velocity  values  at  the  ghost  point  are  determined  by 


U1 

-cos20  +  sin20  -2cos0sin0 

U2 

,vi. 

-2cos0sin0  -sin20  +  cos20 

V2. 

For  the  MacCormack  scheme,  the  normal  component  of  velocity  at 
the  node  point  on  the  surface  is  set  to  zero  and  the 
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Figure  2 .  Body  Surface  "No-Flow  Through"  Boundary  Condition 

tangential  component  is  set  equal  to  the  value  at  the  first 
node  point  off  the  surface. 


To  obtain  the  pressure  at  the  body  surface,  the  normal 
momentum  equation  is  used  as  recommended  by  Pulliam  and 
Steger13: 

-pC7c(tixu£+tiyv{)  =  (nx$x+Eyny)P;  +  (Tix+fly)Pn 

=  Pnh*  +  Vy 


where  p  is  the  pressure  normal  to  the  body  surface.  Second- 


order  accurate  central-differencing  is  used  for 
both  $  and  r\  direction  derivatives  and  Uc  ,the  contravariant 
velocity,  is  given  by 

Uc  =  ixu  +  tyv  (17) 

An  adiabatic  wall  condition  and  the  equation  of  state  are 
used  to  calculate  the  density.  Even  though  the  specification 
of  the  temperature  gradient  at  the  body  surface  is 
inconsistent  wich  the  Euler  equations,  Pulliam  and  Steger13 
report  that  it  has  been  used  by  others  with  no  apparent 
degradation  of  the  solution.  A  recent  paper  by  Driver  and 
Beran14,  using  a  TVD  scheme  similar  to  the  one  developed  in 
this  study,  used  this  same  boundary  condition  for  the  blade 
surface  in  a  high-work,  low  aspect  ratio  turbine  flow 
numerical  simulation.  Their  surface  results  were  all  well 
within  a  few  percent  of  experimental  data.  These  results 
would  seem  to  justify  the  use  of  this  particular  boundary 
condition.  The  total  energy,  Et,  is  calculated  from  the 
surface  values  of  density,  pressure,  and  velocity. 

For  the  TVD  scheme,  the  normal  momentum  equation  is 
applied  at  the  first  cell  center  off  the  body  surface  to 
obtain  values  at  the  ghost  point.  For  the  MacCormack  scheme, 
the  equation  is  applied  at  the  first  node  point  off  the 
surface  to  obtain  values  at  the  surface. 

A  reflection  method  was  also  investigated  for  the  TVD 
scheme  whereby  density  and  total  energy  at  the  ghost  points 
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were  set  equal  to  the  values  at  the  first  cell  center  off  the 
body  surface.  This  technique  slightly  overpredicted  pressure 
at  the  surface  compared  to  pressure  predicted  by  the 
application  of  the  normal  momentum  equation.  Furthermore, 
Roache15  has  shown  that  the  reflection  method  is  an 
inconsistent  boundary  condition  for  a  curved  surface  and 
should  only  be  used  for  flat  surfaces.  Anderson  et  al8  also 
stipulate  that  the  reflection  boundary  condition  "is  very 
inaccurate  on  bodies  in  regions  with  high  curvature." 

2 . 4  Initial  Conditions 

Many  numerical  schemes  are  extremely  sensitive  to  initial 
conditions  and  their  effects  on  the  transient  or  state-state 
solution  is  an  area  of  ongoing  research.  Care  must  be  taken 
to  apply  reasonable,  physically  realistic  initial  conditions 
to  start  calculations16.  For  these  reasons,  detailed 
discussions  of  the  initial  conditions  for  the  steady-state  and 
the  unsteady  computations  will  be  presented  in  the  following 
sections . 

2.4.1.  Steady-State,  Blunt-Bodv  Problem.  For  numerical 
computations  of  the  blunt-body  problem  in  hypersonic  flows, 
many  numerical  procedures  use  theory  to  predict  the  shock 
standoff  distance  and  bow-shock  shape,  and  then  apply  the 
Rankine-Hugoniot  relations  and  modified  Newtonian  theory  to 
provide  a  smooth  initial  state  of  flowfield  quantities.  In 


18 


this  study,  a  first-order  TVD  scheme  was  employed  to  establish 
an  initial  condition  for  the  second-order  TVD  and  MacCormack 
schemes  for  a  number  of  reasons.  First  of  all,  as  part  of  the 
evolutionary  development  of  the  second-order  TVD  scheme,  a 
first-order  TVD  scheme  was  constructed.  It  was  a  simple 
procedure  to  incorporate  that  solution  as  an  initial  state  for 
the  second-order  schemes.  Secondly,  that  initial  state 
contained  a  discontinuity  in  pressure  at  the  stagnation  point 
and  an  unsmooth  flowfield  as  seen  in  the  Mach  contour  values 
(which  will  be  shown  in  Section  V) .  This  provided  a  good  test 
of  the  robustness  of  the  second-order  TVD  and  MacCormack 
schemes.  Finally,  theoretical  solutions  for  an  initial  state 
of  a  general  hypersonic  blunt-body  problem  may  not  always  be 
readily  available  to  a  design  engineer.  The  first-order  TVD 
scheme  is  very  robust  and  can  provide  a  rough  prediction  of  a 
blunt-body  flowfield.  Of  course,  any  numerical  solution 
should  be  compared  to  related  theoretical  predictions  to 
ensure  that  qualitative  features  are  correctly  modelled. 

Within  the  first-order  TVD  scheme,  the  flowfield  is 
initialized  to  freestream  conditions  throughout  and  then  run 
for  500  iterations  (time  steps)  with  a  CFL  number  of  0.9. 
This  technique  is  sometimes  referred  to  as  an  impulsive  start. 
Peyret  and  Taylor16  report  that  Grossman  and  Moretti  found  that 
an  impulsive  start  tended  to  make  transonic  calculations 
unstable  while  a  gradual  or  ramped  start  worked  reasonably 
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well.  This  technique  has  also  been  employed  in  Program  EAGLE17 
for  transonic  calculations  in  the  form  of  gradual  application 
of  the  "no-flow  through"  body  surface  boundary  condition.  For 
the  supersonic/hypersonic  numerical  simulations  in  this  study, 
the  impulsive  start  of  the  first-order  TVD  scheme  was  very 
stable.  The  impulsive  start  procedure  allowed  a  rough 
development  of  the  salient  flow  features  -  such  as  shock 
standoff  distance  and  bow-shock  shape  -  to  take  place.  As 
mentioned  above,  this  solution  constituted  the  initial  state 
used  by  the  second-order  TVD  and  MacCormack  schemes.  The 
second-order  TVD  scheme  was  normally  run  for  500  to  1000 
iterations  with  a  CFL  number  of  0.5  before  the  convergence 
criteria,  to  be  discussed  in  the  next  section,  were  satisfied. 
For  the  MacCormack  scheme,  at  Mach  numbers  greater  than  3.0, 
it  was  necessary  for  stability  to  keep  the  CFL  number  at  or 
below  0.1.  In  addition,  the  smoothing  coefficients  for  the 
MacCormack  scheme  were  set  to  4.0.  Comparison  plots  of  the 
initial  state  from  the  first-order  TVD  scheme  and  final 
converged  solutions  from  the  second-order  TVD  and  MacCormack 
schemes  will  be  presented  in  Section  V. 

2.4.2.  Unsteady  Shock  Impingement  Problem .  For  the 
numerical  simulation  of  the  unsteady  shock  impingement  on  a 
sphere,  an  initial  shock  location  was  assumed  at  x  =  -0.75  for 
the  Mach  2.89  case.  Flowfield  values  are  initialized  to  pre- 
and  post-shock  conditions  and  the  unsteady  solutions  are 
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impulsively  started  from  that  state.  Because  the  dynamic 
process  of  the  shock  striking  and  passing  over  and  beyond  the 
sphere  occurs  on  the  order  of  milliseconds,  the  incident  shock 
wave  can  be  approximated  with  good  accuracy  by  a  normal  shock 
wave18.  As  mentioned  in  Section  2.3.2,  it  was  also  necessary 
to  track  the  shock  along  the  outer  boundary  and  impose  pre- 
and  post-shock  conditions  there  as  the  shock  propagated  in 
time . 

2 . 5  Theoretical  Predictions 

To  ensure  that  the  numerical  simulations  accurately  model 
the  qualitative  features  of  the  steady-state,  blunt-body 
problem  and  the  problem  of  unsteady  shock  impingement, 
comparisons  will  be  made  to  theoretical  predictions  in  the 
results  sections.  The  next  few  sections  present  the  theory. 

2.5.1.  Shock  Standoff  Distance.  The  inviscid  flowfield 
in  the  stagnation  region  between  the  shock  and  the  body 
surface  for  steady,  axisymmetric  flow  about  a  blunt-body  at 
hypersonic  Mach  numbers  is  governed  by  the  hypersonic  boundary 
layer  equations.  Furthermore,  in  the  stagnation  region,  the 
detached  bow-shock  wave  is  nearly  normal  and  hence,  the  Mach 
number  is  subsonic  and  the  flow  is  considered  incompressible. 
Shock  standoff  distance  formulas  are  derived  as  functions  of 
the  density  ratio  and  the  radius  of  the  shock.  The  Rankine- 
Hugoniot  relation  for  the  density  ratio  e  =  p„/ps  across  the 
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normal  shock,  in  the  limiting  case  of  infinite  Mach  number  and 
for  a  value  of  the  ratio  of  specific  heats,  y  =  1.4,  becomes 


€  =  3-^-  =  0.167  (18) 

y+1 

As  reported  by  Reshotko19,  Hayes  derives  the  shock  standoff 
distance,  50  ,  under  these  assumptions  as 


e 

1  +  y/2e 


(19) 


where  Rs  -  radius  of  the  shock.  Additionally,  Truitt20  derives 
the  following  formula  for  axisymmetric  blunt-body  flow: 


-  ci  -i  Wi  -  d-e)2 
Rs  (1-e)2 


(20) 


It  is  readily  apparent  that  an  increase  in  density  in  the 
stagnation  region  will  decrease  the  shock  standoff  distance 
for  both  theoretical  formulas.  This  behavior  will  be  shown  in 
the  tabulated  results  in  Section  V.  For  the  Mach  numbers 
considered,  the  density  ratio  was  calculated  from  the  normal 
shock  relations  as  tabulated  in  NACA  Report  113521.  Anderson 
has  an  excellent  discussion  and  derivation  of  the  limiting 
form  of  the  hypersonic  shock  relations  in  Reference  10. 
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2.5.2.  Modified  Newtonian  Theory.  Also  in  Reference  10, 


Anderson  provides  the  following  historical  perspective: 

In  Propositions  34  and  35  of  his  Principia . 
first  published  in  1687,  Newton  modeled  a  fluid 
flow  as  a  stream  of  particles  in  rectilinear 
motion,  much  like  a  stream  of  pellets  from  a 
shotgun  blast  which,  when  striking  a  surface, 
would  lose  all  their  momentum  normal  to  the 
surface  but  would  move  tangentially  to  the 
surface  without  loss  of  tangential  momentum. 

This  led  to  the  famous  newtonian  sine-squared  law  for 

pressure  coefficient 

Cp  =  2sin20  (21) 

where  0  is  the  local  deflection  angle  of  the  surface. 
Anderson  also  reports  the  Lester  Lees  modification  to 
Newtonian  thee  y  which  gives  the  pressure  coefficient  as 


C  =  C  sin20 

P  P  MAX 


(22) 


Cc  is  the  maximum  value  of  the  pressure  coefficient, 
evaluated  at  a  stagnation  point  behind  a  normal  shock  wave 
given  as 


'  Pmax 


(23) 


where  p,  is  the  total  pressure  behind  the  normal  shock  wave 

''stag 

at  the  freestream  Mach  number.  For  the  Mach  numbers 
considered,  pt  was  determined  from  the  tabulated  normal 

u stag 

shock  relations  in  NACA  Report  113521  . 


23 


2.5.3.  Entropy  Jump  Condition.  The  theory  presented  in 
this  section  is  from  Reference  7.  From  the  second-law  of 
thermodynamics,  the  entropy  always  increases  across  a  shock 
wave.  With  the  assumption  of  a  calorically  perfect  gas,  in 
which  the  coefficients  of  specific  heat  are  constant,  the 
first-law  of  thermodynamics  provides  the  increase  in  entropy 
as 


T  P 

As  =  cDln(— ^22)  -  Rln(— 


(24) 


where  As  =  the  change  in  entropy  from  freestream  to  that 
behind  the  shock,  Tt  -  total  temperature,  R  =  gas  constant, 
Pc  =  total  pressure  and  the  subscripts  Stag  and  °°  represent 
stagnation  and  freestream  quantities,  respectively. 

For  a  stationary  normal  shock  and  a  non-chemically 
reacting  gas 


and  Eq  (24)  becomes 


Tr  =  Tm 

zStag 


„  m( -Eh.) 


(25) 


R 


Pc 


(26) 


Stag 


As  mentioned  in  the  previous  section,  NACA  Report  113521  values 
were  used  to  determine  the  pressures  and  hence  the  theoretical 
entropy  jump  condition. 
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2.5.4.  Shock  Impingement  on  a  Sphere.  Exact  theoretical 
solutions  of  the  dynamics  of  a  planar  moving  shock  wave 
striking  and  passing  over  a  sphere  do  not  exist.  Approximate 
theories  due  to  Whitham  and  to  von  Neumann  are  presented  by 
Bryson  and  Gross22  and  Heilig23,  respectively.  Their 
approximate  natures  only  reveal  the  gross  structure  of  the 
shock  diffraction  process.  Thus,  only  qualitative  features 
from  these  theories  can  be  described.  The  following 
qualitative  description  of  the  time  evolution  of  the  shock- 
diffraction  process  for  a  cylinder  or  a  sphere  is  taken  from 
Young  and  Yee18,  Bryson  and  Gross22,  Heilig23,  Yee  and  Kutler24, 
and  Takayama25. 

The  shock  structure  at  two  instances  of  time,  t,  and  t2, 
after  initially  striking  the  sphere,  is  shown  in  Figure  3. 
This  illustration  is  for  an  incident  shock  of  low  Mach  number 
with  single  Mach  reflection.  When  the  planar  moving  shock 
hits  the  sphere,  regular  reflection  occurs  at  the  shock 
impingement  point.  The  reflected  shock  emerges  as  a  highly 
curved  bow-shock.  The  shock  curvature  induces  an  entropy 
gradient  downstream  which,  in  turn,  introduces  vorticity  into 
the  initially  irrotational  flowfield.  As  the  shock 
impingement  point  propagates  around  the  body,  the  reflection 
process  transitions  from  regular  reflection  to  Mach 
reflection.  Depending  on  the  initial  strength  of  the  shock 
wave,  complex  and  double  Mach  reflection  shock  structures  are 
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Figure  3.  Shock  Diffraction  Process  at  Two  Instances  in 
Time 


also  possible  during  the  transition  process.  For  single  Mach 
reflection,  a  triple  point  forms  and  the  incident  shock  no 
longer  touches  the  body.  Three  waves  emanate  from  the  triple 
point: 

1)  a  Mach  stem  (strong  shock)  which  impinges  normal  to 
the  body  surface, 

2)  a  slip  surface  or  shear  layer  which  strikes  the  body 
and  results  in  a  vortical  singularity  (nodal  point 
of  streamlines) ,  and 

3)  the  reflected  shock  which  moves  away  from  the 
sphere. 

Additionally,  a  stagnation  point  (saddle  point  of  streamlines) 
exists  at  the  lines  of  axisymmetry,  both  fore  and  aft  of  the 
sphere . 

For  double  Mach  reflection,  an  irregular  reflection 
occurs  in  the  vicinity  of  the  triple  point  which  causes  a 


26 


Figure  4.  Double  Mach  Reflection,  Mach  =  3.36 


splitting  into  two  triple  points  -  usually  for  incident  shocks 
greater  than  about  Mach  1.4.  Figure  4  presents  a  Schlieren 
photograph  from  Reference  23  that  illustrates  this  phenomena. 
Whereas  the  Mach  stem  r  hock  is  newly  created  and  was  not 
present  before,  the  reflected  shock  of  the  irregular 
reflection  phase  will  interact  with  the  remaining  reflected 
shock  of  the  regular  reflection  phase.  It  cannot  be  assumed 
that  both  shocks  will  go  together  either  in  their  inclinations 
or  in  their  strengths.  Hence,  their  intersection  appears  as 
a  small  kink  which  starts  at  a  new  triple  point  and  ends  at 


the  original  triple  point.  If  two  shocks  with  unequal  slopes 
interact  in  two  dimensions,  a  new  shock,  "the  second  Mach 
shock,"  and  a  new  slip  line,  "the  second  slip  line",  are 
created.  The  high-resolution  capability  of  the  TVD  scheme 
captures  this  phenomena  as  will  be  shown  in  Section  VI. 

2 . 6  Convergence  Criteria 

One  of  the  convergence  criteria  for  the  TVD  and  the 
MacCormack  scheme  was  the  Euclidean  L2  norm  given  as 

L2  =  (££  (U?;j  -  Ufj)  2)  ^  (27) 

J-l  -2-1 

where  I  and  J  correspond  to  the  maximum  number  of  grid  lines 
in  the  circumferential  and  radial  directions,  respectively, 
and  where  Ui  j  is  the  vector  of  unknowns  given  by  Eq  (2)  .  For 
the  steady-state  computations,  a  reduction  of  the  L2  norm  by 
three  to  four  orders  of  magnitude  from  that  at  the  initial 
state  was  considered  a  converged  solution.  Additionally,  the 
variation  of  the  stagnation  pressure  and  the  shock  standoff 
distance  was  monitored.  These  values  stabilized  quickly,  and 
the  L2  norm  was  the  final  convergence  criterion.  For  the 
unsteady  shock  impingement  computations,  the  L2  norm  was 
monitored  and  stable  values  were  of  order  one. 
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Ill .  Experimental  Data 


For  the  steady-state  blunt-body  problem,  the  experimental 
data  consists  of  graphical  results  for  shock  standoff  distance 
and  tabulated  results  for  wave  shape  and  surface  pressure 
distributions  about  the  nose  of  AGARD  Model  E,  a  hemisphere 
cylinder  configuration.  The  tests  were  conducted  in  the  von 
Karman  Gas  Dynamics  Facility  (VKF)  at  the  Arnold  Engineering 
Development  Center  at  Mach  numbers  3.0,  4.03,  5.06,  6.03,  and 
8 . 1  at  zero  angle  of  attack.  The  experimental  shock  standoff 
distance  and  wave  shape  coordinates  were  obtained  in  the 
report  from  scaled  Schlieren  photographs.  Model  pressures 
were  measured  with  a  system  of  15- ,  5-,  and  1-psid  transducers 
connected  to  orifices  along  the  body.  Figure  5  shows 
excellent  agreement  between  the  experimental  surface  pressures 
and  that  predicted  by  modified  Newtonian  theory,  Eq  (22) . 
Complete  details  of  the  test  can  be  obtained  in  Reference  26. 

For  the  unsteady  shock  impingement  on  a  sphere,  the 
experimental  data  for  visual  comparison  consists  of  Schlieren 
photographs  of  the  diffraction  patterns  at  a  time  when  the 
incident  shock  had  passed  beyond  the  aft  end  of  the  sphere. 
"An  English  table-tennis  ball  of  1  in.  diameter  filled  with 
Wood's  metal  and  suspended  in  the  test  section  from  8  nylon 
strings  was  used  for  a  sphere22."  The  tests  were  conducted  at 
Harvard  University  with  a  shock  Mach  number  of  2.89 
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Figure  5.  Surface  Pressure  Distributions,  Modified 

Newtonian  Theory  and  AEDC  Data,  Mach  =  3.0 


in  a  shock  tube  in  air.  General  features  of  the  shock 
diffraction  process,  as  discussed  in  Section  2.5.4,  are 
readily  discernible  in  the  Schlieren  photographs  to  be  shown 
in  Section  VI. 
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IV.  Numerical  Algorithm  Development 


An  explicit,  second-order  accurate,  upwind,  TVD  scheme 
originally  developed  by  Harten  and  modified  by  Yee1  for  two- 
dimensional  flow  was  extended  to  axisymmetric  flow  in  this 
study.  An  excellent  review  of  TVD  schemes  and  their 
development  can  be  found  in  Reference  1.  Strang  time¬ 
splitting  was  used  to  maintain  second-order  time  accuracy  in 
which  the  local  characteristic  solution  t/mi  is  computed 
from  Un  by 


■♦■1  r  to/ 2  r  to  j-  to/ 2  rj  n 

Ui.j  =  Lr\  Ui,j 


where  £  =  i A£,  ti  =  j'Ar|,  h  =  At  and 


(28) 
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-  E  i  .) 


/  rn* 

(F.  .  x 


-  F 


At 

At]  £'j 


(29) 


This  sequence  of  operators  is  consistent  if  the  sum  of 
the  time  steps  for  each  of  the  operators  are  equal.  Second- 
order  time  accuracy  is  obtained  if  the  operators  are  applied 
in  a  symmetric  sequence8.  Furthermore,  the  time-splitting 
allows  one  to  handle  the  homogeneous  part  of  the  governing 
partial  differential  equations  and  the  source  term  in  separate 
steps1.  The  functions  E .  _i  .  and  F.  .  j.  are  the  numerical 
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fluxes  in  the  £  and  i)  directions,  respectively.  They  are 
evaluated  at  cell  interfaces  and  can  be  written  as 


Ei  +  1  ~  \  t&i  +  &ui  + 

2  ^  2 


(30) 


where  (dropping  subscripts  for  notational  convenience) 


#  =  y„£  -  xnF 

&=  -y^E  +  x5F 


(31) 


This  formulation  has  the  property  of  maintaining  the 
freestream,  thereby  eliminating  the  problem  of  delineating 
physical  flowfield  structure  from  aberrations  due  to  the 
grid27.  In  the  discussion  of  boundary  conditions  in  Section 
2.3.1,  it  was  stated  that  the  flux  at  the  line  of  axisymmetry 
was  set  to  zero.  For  the  front  stagnation  line  with  i  =  2  (i 
=  1  is  the  ghost  point)  the  £  direction  numerical  flux  becomes 


Ei  =  1  (32) 

2  z  2 

with  no  change  in  the  i?_3  numerical  flux.  In  the  unsteady 

2 

computations,  for  the  rear  stagnation  line  with  i  =  1-1  (i  = 
I  is  the  ghost  point)  the  $  direction  numerical  flux  becomes 


-|(*{*) 


(33) 


with  no  change  in  the  Ex  ^_  numerical  flux.  Roe's  averaging 

2 

for  a  perfect  gas11  is  used  to  assess  the  eigenvalues  and 
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eigenvectors  at  cell  interfaces,  because  it  has  the 
computational  advantage  of  resolving  stationary 
discontinuities28.  The  Roe  averaged  dependent  variables  are 
given  by 


Du  ■  , 

i 

+ 

£>  + 
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+ 

D  + 
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+ 

-  <v-i>  '"i.i  -  j 

2  2*2 
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)] 


where 


D  ^P  i*\/ Pi 


(34) 


(35) 


The  vector  elements  of  for  the  TVD  scheme  are 


+  (9i\ i  +  9i)  ~  +  2  (36) 

2  2  2  2  2 

and  I  =  1  to  4  is  consistent  with  Eq  (2)  and  where,  with 
A=At/Ax 


i  * 


_1 

2 


(37) 
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ai  +  j_  is  the  difference  of  the  characteristic  variables 

2 


in  the  $  direction, 


«i.J:  =  R-UiU^  -  Ui] 

o  -*■  "r  T" 


(38) 


The  components  of  a  are  given  by 


a 
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“U 
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*4 


(aa  -  bb)  /2 
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(aa  +  bb)  /2 
cc 
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where 


aa 


bb 


=  [A  .  i  e  + 


1  2 


+  V 


1+  2 


lkhAi.±mx  ~  (khlUu±  +  +  kh2*i.±Wy 


CC  =  knAi>±my  +  ^kh2^i.l  -  -  kh2^Ul^X 


and 


Ai.Az  =  zi*i  "  zi 
2 


(40) 


The  difference  of  the  characteristic  variables  in  the  rj 
direction  is  constructed  similarly. 
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The  function  y1  x  is 

i*-=- 


where 


and 
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|z|  ^  8f 

|z|  <  8f 


(41) 


(42) 


(43) 


The  functions  o  and  are  referred  to  as  the  effective 
numerical  viscosity  coefficient  and  the  numerical  viscosity 
coefficient,  respectively4.  The  functional  form  of  i|r  was 
developed  by  Harten  in  Reference  4.  It  serves  as  an  entropy 
correction  to  \z\  near  0  to  prevent  an  entropy  violation  which 
occurs  when  Eq  (42)  vanishes.  One  can  view  the  size  of  the 
entropy  correction  function,  8f  ,  as  a  measure  of  the  amount 
of  numerical  dissipation  added.  be  =  0  is  the  least 
dissipative,  and  the  larger  the  bf  the  more  dissipative  the 
scheme  becomes.  The  numerical  simulation  of  hypersonic 
blunt-body  flows  is  especially  sensitive  to  the  magnitude  of 
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the  entropy  correction  function.  As  reported  in  Reference  1, 
Yee  suggests  the  following  form  employed  in  this  study: 


(6  f) 


=  »p(|(uc) 


'4' 


<4 


cu±) 

2 


(44) 


where  5p  is  the  entropy  correction  parameter  and  uc  and  vc  are 
the  contravariant  velocities.  Additionally,  half  of  the  sound 
speed  is  from  the  £  direction  and  the  other  half  is  from  the 
r;  direction.  The  function  5f  depends  on  the  spectral  radius 
of  the  Jacobian  matrices  of  the  fluxes  and  is  very  useful  in 
terms  of  stability  and  convergence  rate.  Yee  recommends 
0.05  ^  5p  ^  0.25  for  4  s  Mm  <,  25 .  For  the  nondimensional 
quantities  employed  in  this  study,  a  consistent  value 
for  8p  of  5.0  was  used  unless  indicated  otherwise.  It  should 
be  pointed  out  that  smaller  values  slow  down  convergence  for 
the  steady-state  problem.  As  mentioned  previously,  larger 
values  of  8f  add  more  numerical  dissipation.  An  investigation 
of  these  effects  was  conducted  and  3.5  £  8p  s  5.5  was  found  to 
work  quite  well.  Values  lower  than  5.0  dramatically  slowed 
down  the  convergence  rate  with  only  a  slight  increase  in  the 
resolution  of  the  shock  as  will  be  shown  in  Section  V.  For 
instance,  after  1000  iterations,  a  value  of  5.0  reduced 
the  L2  norm  three  orders  of  magnitude.  A  value  of  4.0  for  the 
same  number  of  iterations  did  not  reduce  the  L2  norm  a  single 
order  of  magnitude. 


36 


The  function  gf  is  the  'limiter'  function  which  also 
controls  the  amount  of  numerical  dissipation  added. 
Hypersonic  blunt-body  flows  with  strong  shock  waves  require 
the  proper  selection  of  limiters.  The  limiter  function  has 
been  expressed  in  a  number  of  different  ways1.  This  study 
incorporates  limiters  appropriate  for  the  characteristic 
fields  under  investigation  as  recommended  by  Yee1 .  The 
nonlinear  fields  correspond  to  1  =  1  and  1=3  and  these  waves 
are  either  shocks  or  rarefaction  waves.  For  the  nonlinear 
fields,  a1R1  *  0,  use 

9t  =  (a1  1a1  1  +  |aJ  1a1  i|)/(aJ  a  +  a1  1 )  (45) 

i  +  _  1  -  —  l  *■  — i-  7  -  —  14-  -t-  1-  \  / 

2  2  2  2  2  2 

The  linearly  degenerate  fields  correspond 
to  1  =  2  and  1=4  and  are  uniquely  contact  discontinuities6. 
For  the  linearly  degenerate  fields,  a1R1  =  0,  use 

9i  =  minmod{a1._1 ,  )  (46) 

where  the  minmod  function  of  a  list  of  arguments  is  equal  to 
the  smallest  number  in  absolute  value  if  the  list  of  arguments 
is  of  the  same  sign,  or  is  equal  to  zero  if  any  arguments  are 
of  opposite  sign.  Limiters  other  than  Eq  (46)  can  produce 
sharper  discontinuities,  but  are  not  as  robust.  For  the 


37 


unsteady  shock  impingement  problem,  the  following  limiter  was 
used  for  the  linearly  degenerate  fields 

g/=S  •  max  [0,min  (2  |ai  ^  \ ,  S  •  a 1  ^ )  ,  min  (|a^  +  j_|,25-ai_i)]  (47) 

2  '2  +  2  2 

where 

S  =  sgn (a1  x)  (48) 

1  +  -=-  V  7 

A  numerical  investigation  of  the  differences  in  resolution  of 
flowfield  discontinuities  with  these  limiters  will  be 
presented  in  Section  VI.  Finally,  note  that  zeroth-order 
extrapolation  is  used  to  obtain  limiter  values  at  ghost 
points . 
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V.  Steady-State  Solutions  of  the  Euler  Equations 


Steady-state  solutions  obtained  with  the  TVD  and  the 
MacCormack  schemes  are  presented  in  this  section. 

5 . 1  Computational  Grid 

Steady-state  solutions  were  computed  using  the  grid  shown 
in  Figure  6.  The  grid  consisted  of  51  points  in  the 
circumferential  direction  (I  -  51)  and  25  points  in  the  radial 
direction  (J  =  25)  .  The  body  radius,  RB,  was  scaled  to  a 
nondimensional  unit  of  l.  The  spacing  in  the  circumferential 
direction  at  the  line  of  axisymmetry  was  .02  RB  ,  and  the 
spacing  in  the  radial  direction  was  .01  RB.  A  geometric 
progression  was  used  to  control  the  radial  spacing  from  the 
body  surface  to  the  outer  boundary  and  for  the  circumferential 
spacing  from  the  line  of  axisymmetry  to  the  rear  outflow 
boundary.  This  allowed  for  clustering  of  node  points  in  the 
stagnation  region.  The  location  of  the  far-field  boundary  was 
set  to  -1.4  Rg  along  the  line  of  axisymmetry  and  an  outer 
boundary  with  a  parabolic  shape  was  used. 

5 . 2  Shock  Standoff  Distance  and  Wave  Shape 

To  evaluate  the  performance  of  the  TVD  and  MacCormack 
schemes,  comparisons  of  shock  standoff  distance  and  bow-shock 
shape  were  made  to  theory  and  to  experimental  data.  To 
determine  the  shock  standoff  distance  predicted  by  the 
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Figure  6.  Computational  Grid,  51  x  25 


schemes,  plots  of  pressure  along  the  line  of  axisymmetry  were 
used.  An  example  of  such  a  plot  for  the  Mach  3.0  case  is 
shown  in  Figure  7.  Note  that  the  TVD  data  is  from  cell 
centers  and  the  MacCormack  data  is  from  node  points.  These 
expanded  scale  plots  allow  for  an  interpolation  of  the  data  to 


determine  the  shock  standoff  distance.  In  Figure  7,  note  that 
the  TVD  scheme  resolves  the  shock  in  fewer  grid  points  than 
the  MacCormack  scheme.  Typically,  TVD  schemes  do  much  better 
than  this.  It  was  mentioned  in  Section  IV  that  the  entropy 
correction  parameter,  5p  ,  affected  the  amount  of  numerical 
dissipation  added  and,  in  turn,  the  resolution  of  shocks. 
Values  for  6p  less  than  5.0  tended  to  align  the  two  points 
around  the  shock  with  the  pre-  and  post-shock  values.  The  net 
result  was  resolution  of  the  shock  in  a  couple  of  grid  points. 
However,  the  penalty  paid  was  a  dramatic  decrease  in  the 
convergence  rate.  For  this  study,  the  efficient  calculation 
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of  solutions  was  a  primary  concern  to  the  author.  Most  of  the 
numerical  simulations  were  done  interactively  in  a  few 
thousand  iterations  at  most.  Figure  8  shows  a  contour  plot  of 
the  entropy  correction  function  given  by  Eq  (44)  for  the  Mach 
3.0  case.  It  appears  that  the  magnitude  of  the  function  is 
related  to  "position"  in  the  physical  domain.  That  is,  is 
greater  near  the  line  of  axisymmetry  and  decreases  in  the 
vicinity  of  the  shock  where  the  shock  strength  decreases. 
Figure  9  compares  pressure  plots  along  the  line  of  axisymmetry 
from  the  first-order  TVD  scheme  to  the  second-order  scheme  for 
the  Mach  3.0  case.  Note  that  the  first-order  TVD  scheme 
resolves  the  shock  very  well.  However,  the  first-order  scheme 
does  not  do  so  well  in  the  region  between  the  shock  and  the 
stagnation  point. 
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Figure  8.  Entropy  Correction  Function  Contours,  Mach 
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Figure  9.  Pressure  Along  the  Line  of  Axisymmetry ,  First- 
and  Second-Order  TVD  Schemes,  Mach  =  3.0 


Table  1  presents  a  comparison  of  the  calculated  shock 
standoff  distances  to  theory  from  Section  2.5.1  and  AEDC 
experimental  data. 


Table  1 

Shock  Standoff  Distance 

50/*s 


Mach=> 

3 . 0 

4.03 

5.06 

6.03 

TVD 

.  177 

.  144 

.  132 

.  126 

.  151 

.  131 

.  122 

.  117 

1. 

.  155 

.  134 

.  125 

.  119 

AEDC 

.  179 

.  138 

.  126 

.  122 

8.1 


.  121 


.  106 


.  107 


The  TVD  and  MacCormack  schemes  compare  quite  well  with 


theory  and 

experimental  data 

over  a  wide 

range 

of 

Mach 

numbers . 

It  should  also  be 

pointed  out 

that 

for 

the 

experimental  data  the  shock  standoff  distances  were  given  in 
graphical  form.  Hence,  the  three  significant  figures  shown 
are  uncertain. 

The  bow-shock  shape  was  also  predicted  quite  well  by  both 
the  TVD  and  MacCormack  schemes.  Figures  10  and  11  present 
Mach  contours  -  with  the  shock  coordinates  from  the  tabulated 
experimental  data  superimposed  -  from  the  TVD  and  MacCormack 
schemes,  respectively.  The  TVD  scheme  Mach  contours  are  much 
"crisper"  than  those  of  the  MacCormack  scheme,  especially  in 
the  stagnation  region.  Also,  smoother  Mach  contours  at  the 
outflow  boundary  are  obtained  with  the  TVD  scheme  as  shown  in 
Figure  10.  This  graphically  illustrates  the  point  made  in 
Section  2.3.3  concerning  two-point  extrapolation  as  opposed  to 
a  "no-change"  condition  as  the  outflow  boundary  condition. 

Figure  12  presents  Mach  contours  from  the  TVD  solution  at 
a  Mach  number  of  8.1.  Crmpare  this  to  Figure  10  and  the 
statement  from  Anderson10: 

...as  Mm  increases,  the  shock  wave  moves  closer  to  the 
body  and  the  sonic  points  on  both  the  shock  and  the  body 
move  closer  to  the  centerline  -  all  standard  physical 
behavior  for  blunt-body  flows.  Furthermore,  observe 
that,  as  M„  increases,  the  sonic  point  on  the  shock  moves 
down  faster  than  the  sonic  point  on  the  body,  and  thus 
the  sonic  line  actually  rotates  in  a  counterclockwise 
fashion  as  the  Mach  number  increases. 


4  5 
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Figure  11.  Mach  Contours, 


MacCormack  Scheme,  Mach  = 


4.03 
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Figure  12.  Mach  Contours,  TVD  Scheme,  Mach  =8.1 

Furthermore,  for  two-dimensional  flow  about  a  cylinder,  the 
angle  that  the  sonic  line  makes  with  the  body  surface  remains 
acute  -  no  matter  how  high  the  Mach  number.  For  axisymmetric 
flow  about  a  sphere,  that  angle  is  acute  at  low  Mach  numbers 
and  becomes  obtuse  for  Mach  numbers  approximately  greater  than 


3.0.  Figures  10  and  12  demonstrate  that  the  axisymmetric  TVD 
scheme  correctly  models  this  physical  behavior. 

5 . 3  Pressure  Distribution 


To  assess  the  performance  of  the  TVD  and  MacCormack 
schemes,  numerically  predicted  surface  pressure  distributions 
are  compared  to  experimental  data. 


Figure  13  presents  pressure  distribution  calculations 
from  the  first-  and  second-order  TVD  schemes  compared  to  the 
experimental  data  at  Mach  3.0.  The  solution  procedure  was; 
run  the  first-order  TVD  scheme  for  500  iterations  at  a  CFL 
number  of  0.9,  use  that  solution  as  an  initial  state  for  the 
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second-order  TVD  scheme,  and  run  the  second-order  TVD  scheme 
for  as  many  iterations  necessary  (usually  500  to  1000)  at  a 
CFL  number  of  0.5  until  convergence.  Convergence  criteria 
were  presented  in  Section  2.6.  In  Figure  13,  the  first-order 
TVD  scheme  solution  has  a  pressure  oscillation  in  the 
stagnation  region  which  underpredicts  the  stagnation  pressure 
and  the  scheme  overpredicts  the  pressure  distribution.  The 
second-order  TVD  scheme  solution  completely  damps  out  the 
oscillatory  behavior  and  there  is  excellent  agreement  with  the 


experimental  data.  Similar  results  were  obtained  for  Mach 


numbers  of  4.03,  5.06,  6.03,  and  8.1.  Figure  14  presents  a 
pressure  distribution  calculation  from  the  MacCormack  scheme 
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compared  to  the  experimental  data  at  Mach  3.0.  This  solution 
was  obtained  in  10,000  iterations  with  a  CFL  number  of  0.5. 
It  was  not  possible  to  reduce  the  L2  norm  convergence 
criterion  as  low  as  that  obtained  with  the  TVD  scheme.  At 
Mach  numbers  greater  than  3.0,  it  was  necessary  to  reduce  the 
CFL  number  to  no  more  than  0.1  for  stability.  Furthermore, 
after  a  few  thousand  iterations,  an  unphysical  recirculation 
region  would  develop  in  the  stagnation  region,  resulting  in  an 
accentuated  bulge  in  the  bow-shock  shape  and  large  pressure 
oscillations  about  the  stagnation  point.  This  type  of 
behavior  in  the  MacCormack  scheme  for  the  blunt-body  problem 
has  been  reported  by  Shang  and  Josyula6.  In  fact,  the 
accentuated  bulge  looks  like  an  unsteady  shock  boundary  layer 
interaction  ahead  of  a  forward  facing  step  as  reported  by 
Saida  et  al29.  The  causes  for  this  behavior  are  not  completely 
understood.  To  eliminate  the  possibility  that  the  initial 
state  caused  this  behavior  in  the  MacCormack  scheme,  a  smooth 
initial  state  was  developed  as  described  in  Section  2.4.1. 
Unfortunately,  the  recirculation  region  also  developed  from 
this  initial  condition.  An  example  of  the  onset  of  this 
behavior  is  shown  in  Figure  15. 
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Second-Order  TVD  Schemes,  Mach  =  4.03 


5 . 4  Entropy  Jump  Condition 

It  is  well  known  that  not  all  algorithms  for  the  Euler 
equations  will  compute  the  physically  correct  solution  and 
that  an  entropy  condition  is  required  to  pick  out  the  relevant 
solution30.  The  level  of  spurious  entropy  produced  in  the 
stagnation  region  of  a  blunt-nosed  configuration  is  a  good 
measure  of  a  numerical  scheme's  accuracy31. 
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Table  2  compares  theoretical  results  to  that  obtained 


with  the  second-order  TVD  and  MacCormack  scheme. 


Table  2 

Entropy  Jump  Condition 
A  s/R 


Mach=> 

3 . 0 

4.03 

5.06 

6.03 

8.1 

Eq (26) 

1.114 

2.001 

2.831 

3 . 539 

4.825 

TVD 

1.095 

1.977 

2 .805 

3.511 

4.791 

%Error 

-1.73 

-1.18 

-  .93 

-  .78 

-  .70 

Mac 

1.133 

2 . 027 

2 . 875 

3.579 

— 

%Error 

+  1.60 

+  1.30 

+  1.55 

+  1.14 

— 

Excellent  agreement  is  obtained  between  the  TVD  scheme 
and  theory.  The  MacCormack  scheme  results  are  also  quite 
good,  but  stagnation  point  pressure  oscillations  were  evident, 
and  the  L2  norm  had  only  been  reduced  two  orders  of  magnitude 
with  the  very  restrictive  CFL  number  necessary  for  stability. 
Further  iterations  resulted  in  diverging  behavior  as  mentioned 
in  the  previous  section. 
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VI .  Dynamic  Solutions  of  the  Euler  Equations 


The  complex  unsteady  flowfield  generated  by  the 
impingement  of  a  planar  shock  on  a  sphere  is  a  good  test  case 
for  assessing  the  time  accuracy  and  shock-capturing  capability 
of  the  second-order  TVD  and  MacCormack  schemes.  A  discussion 
of  the  time  evolution  of  the  shock  diffraction  process  was 
given  in  Section  2.5.4.  Flowfield  patterns,  specifically 
density  and  pressure  contours,  at  different  instances  in  time 
for  the  Mach  2.89  case  will  be  presented  in  the  next  two 
sections.  A  comparison  of  limiter  effects  on  the  resolution 
of  the  flowfield  structure  will  also  be  presented. 
Additionally,  comparison  to  Schlieren  photographs  from  the 
Bryson  and  Gross22  experiments  will  be  shown. 

6 . 1  Computational  Grid 

The  numerical  simulations  were  obtained  on  the 
cylindrical  grid  shown  in  Figure  16.  The  grid  consisted  of 
251  evenly  spaced  points  in  the  circumferential  direction  and 
101  points  in  the  radial  direction.  The  spacing  in  the  radial 
direction  off  the  body  was  .01  Rg  and  a  geometric  progression 
was  used  to  control  the  spacing  from  there  to  the  outer 
boundary.  The  outer  boundary  had  a  radius  of  2  body 
diameters . 
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6 . 2  Pressure  and  Density  Flowfield  Predictions 


The  numerical  simulations  were  started  with  the  planar 
shock  located  at  x  =  -0.75  with  a  shock  Mach  number  of  2.89. 
Pre-  and  post-shock  conditions  are  imposed  to  the  right  and 
left,  respectively,  of  the  shock  location. 

Figures  17  and  18  are  density  and  pressure  contours, 
respectively,  at  selected  times  in  the  shock  diffraction 
process  for  the  TVD  scheme.  The  limiter  used  for  the  linear 
fields  was  given  as  Eq  (47) .  The  linearly  degenerate  fields 
are  uniquely  contact  discontinuities  as  described  in  Section 
IV.  Figures  19  and  20  are  density  and  pressure  contours, 
respectively,  at  approximately  the  same  selected  times  in  the 
shock  diffraction  process,  for  the  MacCormack  scheme. 

Figures  17a,  18a,  19a,  and  20a  show  the  reflected  shock 
propagating  away  from  the  front  of  the  sphere  with  the 
incident  shock  attached  to  the  body  surface  at  the  impingement 
point.  Note  the  smoother  contours  for  the  TVD  scheme  results. 

In  Figures  17b,  18b,  19b,  and  20b,  the  impingement  point 
has  evolved  into  the  split  triple  point  of  a  double  Mach 
reflection.  This  is  barely  discernible  as  the  kink  in  the 
reflected  shock  -  incident  shock  intersection.  Double  Mach 
reflection  involves  double  Mach  stems  and  slip  surfaces 
(contact  discontinuities) .  These  will  be  referred  to  as 
primary  and  secondary  for  the  "upstream"  and  "downstream" 
locations . 
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Figure  17.  Density  Contours,  Second-Order  TVD  Scheme, 
Time  =  a)  .4924,  b)  .6727 ,  c)  .8703 


Figure  17.  Density  Contours,  Second-Order  TVD  Scheme, 
Time  -=  d)  1.0667,  e)  1.2708,  and  f)  1.4765 


Figure  18.  Pressure  Contours,  Second-Order  TVD  Scheme, 
Time  =  a)  .4924,  b)  .6727,  c)  .8703 


Figure  18.  Pressure  Contours,  Second-Order  TVD 
Time  =  d)  1.0667,  e)  1.2708,  and  f) 


Scheme , 
1.4765 
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Figure  19.  Density  Contours,  MacCormack  Scheme 
Time  =  a)  .4901,  b)  .6706,  c)  .8679 


Figure  19.  Density  Contours,  MacCormack  Scheme, 

Time  =  d)  1.0642,  e)  1.2690,  and  f)  1.4747 
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Figure  20.  Pressure  Contours,  MacCormack  Scheme, 
Time  =  a)  .4901,  b)  .6706,  c)  .8679 


ours 

642, 


MacCormack  Scheme, 
e)  1.2690,  and  f)  1.4747 


The  moment  at  which  the  incident  shock  and  primary  Mach 
stem  pass  over  the  top  of  the  sphere  is  a  critical  time  for 
the  numerical  simulations.  The  flow  numerically  overexpands 
and  the  pressure  becomes  negative  in  some  instances.  This 
problem  has  been  reported  for  the  numerical  simulation  of  a 
two-dimensional  flow  about  a  triangular  obstacle  with  a 
rounded  top1  using  the  MacCormack  scheme.  Variations  of  CFL 
number  and  smoothing  coefficients  for  that  case  only  delayed 
the  onset  of  negative  pressure.  In  this  study,  the  MacCormack 
scheme  was  used  with  a  CFL  number  of  0.5  and  values  for  the 
smoothing  coefficients  of  3.0.  Values  of  the  smoothing 
coefficients  of  2.0  did  result  in  the  flow  numerically 
overexpanding  just  aft  of  the  top  of  the  sphere.  The  second- 
order  TVD  scheme  used  approximately  the  same  operating 
parameters  here  as  in  the  steady-state  problem,  i.e.,  CFL 
number  of  0.5  and  5p  =  3.0. 

Figure  17c,  density  contours  from  the  TVD  scheme,  shows 
the  emergence  of  the  primary  slip  surface  just  aft  of  the  top 
of  the  sphere.  In  Figure  17c,  the  slip  surface  has  rolled  up 
into  a  vortex  -  not  to  be  confused  with  the  vortices  that  form 
downstream  of  spheres  due  to  separation.  Figure  19c,  density 
contours  for  the  MacCormack  scheme  at  the  same  time  instant, 
does  not  resolve  the  slip  surface  as  distinctly  as  the  TVD 
scheme.  In  the  pressure  contours  for  both  schemes  of  Figures 
18c  and  20c,  it  is  possible  to  see  where  the  slip  surface 
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impinges  the  body.  At  that  point,  a  local  pressure  minimum 
occurs,  and  the  pressure  contours  encircle  it. 

The  density  contours  in  the  time  intervals  of  Figures  17d 
to  17 f  for  the  TVD  scheme  show  the  emergence  of  the  secondary 
slip  surface.  This  is  discernible  as  a  "sinusoidally  shaped" 
kink  in  the  contours.  This  surface  does  not  roll  up  into 
vortices  like  the  primary  slip  surface  does.  The  density 
contours  of  Figures  19d  to  19f  for  the  MacCormack  scheme  do 
not  resolve  the  secondary  slip  surface. 

Furthermore,  during  this  same  time,  the  secondary  Mach 
stem  rotates  clockwise  down  to  the  body  surface  and  propagates 
upstream  within  the  reflected  shock  -  body  surface  region. 
The  pressure  contours  in  Figures  18d  to  18 f  for  the  TVD  scheme 
and  Figures  20d  to  20f  for  the  MacCormack  scheme  do  not 
resolve  this  phenomena.  Figure  21  is  a  plot  of  the  surface 
pressure  distributions  which  correspond  to  Figures  18d  and  20d 
for  the  TVD  and  MacCormack  schemes,  respectively.  The  TVD 
scheme  shows  a  small  rise  in  the  surface  pressure  at  this 
secondary  Mach  stem,  just  ahead  of  the  shock  impingement 
point,  at  approximately  110  degrees  from  the  front  stagnation 
point.  The  MacCormack  scheme  shows  a  gradual  rise  in  the 
surface  pressure  from  this  point  to  the  shock  impingement 
point.  Note  also  in  Figure  21,  for  the  MacCormack  scheme,  the 
onset  of  pressure  oscillations  at  the  stagnation  point  similar 
to  those  when  the  MacCormack  scheme  was  applied  to  the  steady- 
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state  problem.  These  stagnation  point  oscillations  also 
appear  in  the  TVD  scheme  results  as  the  shock  impingement 
point  moves  further  aft  toward  the  rear  stagnation  point. 
Bennett  et  al32  came  upon  this  same  phenomena  for  shock 
impingement  on  a  cylinder  using  ARC2D  -  an  inviscid,  finite- 
difference  scheme  based  on  the  Beam-Warming  algorithm. 
Bennett  et  al  extended  the  code  to  the  Navier-Stokes  equations 
and  found  that  the  pressure  oscillations  were  removed.  An 
interesting  aside  in  the  Bennett  report  was  that  the 
experimental  data  obtained  by  Pearson  at  the  Ballistic 
Research  Laboratory  for  the  shock  impingement  on  a  cylinder 
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These 


had  "substantial  high  frequency  oscillations." 
oscillations  were  time-averaged  for  the  final  data  set  used  in 
the  comparisons  to  the  inviscid  and  viscous  ARC2D  numerical 
predictions  of  pressure  distributions. 

A  comparison  of  the  effects  of  different  limiters  for  the 
TVD  scheme  can  be  seen  in  the  density  contours  of  Figure  22. 
These  were  obtained  at  approximately  the  same  time  instant  as 
those  of  Figure  18d  using  the  limiter  of  Eq  (46)  for  the 
linearly  degenerate  fields.  This  dissipative  limiter  barely 
resolves  the  primary  slip  surface  and  totally  ignores  the 
secondary  slip  surface.  For  this  reason,  the  limiter  of  Eq 
(47)  is  recommended. 


•2  00  - 1  00  0  00  1  00  2  CO 

X 

Figure  22.  Density  Contours,  Second-Order  TVD  Scheme, 
Limiter  of  Eq  (46),  Time  =  1.1546 

A  comparison  between  the  results  from  the  TVD  scheme  and 
Schlieren  photographs  from  the  Bryson  and  Gross22  experiment 
are  shown  in  the  density  contours  in  Figure  23  and  24.  These 


68 


Figure  23.  Density  Contours,  Second-Order  TVD  Scheme, 
Time  =  1.7510,  ;c  =  .9990 


I 


70 


Figure  24.  Density  Contours,  Second-Order  TVD  Scheme, 
Time  =  2.0527,  x  =  1.3007 


correspond  to  incident  shock  locations  of  approximately  1  and 

1.3  body  diameters  aft  of  the  center  of  the  sphere, 
respectively.  The  eight  radial  lines  in  the  photos  are  nylon 
strings  used  to  support  the  sphere  in  the  test  section.  The 
triple  points,  slip  surfaces,  vortices,  Mach  stems,  and 
reflected  shocks  are  all  favorably  resolved  by  the  TVD  scheme 
when  compared  to  the  Schlieren  photographs. 

6 . 3  Comparison  of  Run  Times  and  Time  Step  Count 

A  comparison  of  the  run  time  for  each  time  step  and  the 
number  of  steps  required  to  reach  the  same  selected  times  for 
the  TVD  and  MacCormack  schemes  will  be  presented  in  this 
section. 

A  Stardent  GS-2000  Graphics  Supercomputer  of  the 
Aeronautics  and  Astronautics  department  of  the  Air  Force 


Institute 

of 

Technology 

was  used 

for 

this 

study. 

The 

numerical 

TVD 

scheme  was 

optimized 

to 

allow 

for  as 

much 

vectorization  as  possible.  The  MacCormack  scheme  has  also 
been  optimized.  For  the  251  x  101  grid,  the  computation  times 
were  approximately  4  seconds  per  time  step  and  2  seconds  per 
time  step  for  the  TVD  and  MacCormack  schemes,  respectively. 
This  •  ,i responds  to  a  data  processing  rate  of  1.58  x  10'4  and 
7.89  x  10  5  seconds  per  grid  point  per  time  step,  respectively. 
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Table  3  compares  the  number  of  time  steps,  time,  and 
percent  time  step  increase  for  the  TVD  scheme  and  MacCormack 
schemes . 


Table  3 

Comparison  of  Time  Step  Count  and  Time 


‘  TVD  Scheme 

MacCormack  Scheme 

%  Count 
Rise 

Count 

Time 

Count 

Time 

200 

.  4924 

258 

.4901 

29.0 

300 

.6727 

386 

.  6706 

28.7 

400 

.8703 

517 

.8679 

29.3 

500 

1.0667 

636 

1.0642 

27.2 

600 

1.2708 

752 

1.2690 

25.3 

700 

1.4765 

869 

1.4747 

24 . 1 

800 

1.6775 

985 

1.6751 

23.1 

For  the  same  CFL  number,  the  MacCormack  scheme  has  a  more 
restrictive  time  step  size  than  the  TVD  scheme. 

Finally,  for  the  251  x  101  grid,  to  reach  a  time  of 
1.6775,  it  takes  approximately  35%  longer  with  the  TVD  scheme 
than  the  MacCormack  scheme.  This  is  still  under  1  hour  and 
quite  reasonable  for  a  numerical  simulation  of  this 
complexity. 
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VII.  Conclusions  and  Recommendations 


An  explicit,  second-order  accurate,  finite-volume,  TVD 
scheme  has  been  developed  for  the  Euler  equations  in 
axisymmetric  form.  A  summary  of  conclusions  from  the  results 
presented  in  this  study  and  recommendations  for  further 
investigation  are  presented  in  the  next  two  sections. 

7 . 1  Conclusions 

Numerical  simulations  of  the  steady-state,  blunt-body 
problem  were  in  excellent  agreement  with  theory  and 
experimental  data  over  a  wide  range  of  Mach  numbers  in  terms 
of  shock  standoff  distance,  bow-shock  shape,  and  surface 
pressure  distributions.  The  calculated  entropy  jump  was  also 
in  excellent  agreement  with  theory  and  experimental  data. 
Comparisons  with  the  MacCormack  scheme  show  that  the  second- 
order  TVD  scheme  better  resolves  the  flowfield  features  and 
has  a  greater  stability  and  robustness  characteristic.  It  was 
shown  that  the  choice  of  the  entropy  correction 
parameter,  6p  ,  affects  convergence  rate  and  resolution  of 
flowfield  features. 

Numerical  simulations  of  the  shock  impingement  problem 
with  both  schemes  show  that  the  TVD  scheme's  ability  to 
resolve  the  complex,  unsteady  interactions  during  a  blast-wave 
impact  on  a  sphere  are  very  good.  Primary  and  secondary 
structures  were  resolved  with  the  proper  limiter  for  the 
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linearly  degenerate  fields.  These  flowfield  structures  were 
not  resolved  as  well  with  the  MacCormack  scheme.  For  the  TVD 
scheme,  comparisons  of  flowfield  structure  to  Schlieren  photos 
form  experiment  were  also  very  good. 

These  results  suggest  that  the  proposed  TVD  scheme  can  be 
a  very  useful  tool  in  the  efficient  calculation  of  accurate 
solutions  for  the  research,  development  and  design  of 
hypersonic  blunt-body  vehicles. 

7 . 2  Recommendations 


Recommendations  for  further  study  of  the  TVD  scheme  that 
are  applicable  to  the  steady-state  blunt-body  problem  and  the 
unsteady  shock  impingement  problem  are: 

1)  a  systematic  study  of  limiters  and  their  effect  on 
resolution  of  discontinuities, 

2)  extension  to  three-dimensional  form  using  the 
axisymmetric  form  as  a  test  bed, 

3)  incorporation  of  equilibrium  and  non-equilibrium 
effects,  and 

4)  extension  to  the  Navier-Stokes  equations. 


Difficulties  anticipated  with  these  recommendations 
include : 

1)  loss  of  time  accuracy  and 

2)  increased  computer  memory  and  CPU  time  requirements 
due  to  increased  number  of  nodes  and  complexity  of 
the  governing  equations. 
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Recommendations  for  further  study  of  the  steady-state 
problem  blunt-body  problem  include: 

1)  automatic  switching  from  a  first-order  to  a 
second-order  scheme  for  an  impulsive  start,  or 

2)  gradual  application  of  boundary  conditions  to  allow 
an  impulsive  start,  and 

3)  further  analysis  of  the  effects  of  the  entropy 
correction  function  and  parameter  on  convergence 
rate  and  resolution. 

Recommendations  for  further  study  of  the  unsteady  shock 
impingement  problem  include: 

1)  simulations  for  different  structures  -  data  from 
nuclear  blast-wave  effects  are  important  for 
national  security  reasons33 

2)  analysis  and  investigation  into  the  pressure 
oscillation  phenomena  (perhaps  the  Navier-Stokes 
equations  will  remove  the  phenomena  as  it  did  for 
Bennett  et  al32)  ,  and 

3)  development  of  better  freestream  boundary  conditions 
to  eliminate  having  to  track  the  shock. 
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Appendix:  Nondimens ional ization 


The  Euler  equations  of  motion  are  nondimensionalized 
following  the  description  by  Beran  in  Reference  34.  The 
freestream  velocity,  Um,  and  the  sphere  radius,  Rg,  are  the 
velocity  and  length  scales,  respectively.  The  freestream 
conditions  are  specified  to  be  scales  for  the  temperature  and 
density.  Consequently, 

u  =  U'U„  V  =  v'U" 

X  =  x'Rb  y  =  y’RB 
t  =  f(Rg/Um) 

T  =  T*T„  p  =  P*P„ 

P  =  P*p«y«  e  =  e*ul 


The  Mach  number,  =  Ujc„,  where  cm  is  the  freestream  speed 
of  sound,  is  used  in  the  boundary  conditions.  A  more 
convenient  form  of  the  Mach  number  is 


Ml 


Y 


where 


qi  =  ut 


+  v_ 
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The  equation  of  state  (perfect  gas)  ,  p  =  p RT,  is  written  in 
non-dimensional  form  as 


P’p.C/i  =  p'p  JiT’Tm 
p'ui  =  p  'T'(RTJ 


P*  = 


p'T’ 
Y  Ml 


UZ 

PT( - -) 

yMZ 


Note  that  in  the  freestream 


p*  =  (Y Ml)  - 1 


The  pressure  is  eliminated  in  the  equations  of  motion  through 
another  approach: 


p  =  9r~~-  ~  p  (y-D  e 
p*  =  (Y-l)p’e* 


Note  that  in  the  freestream, 

e*  =  (  (y-1)  Y Mi) 
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The  definition  of  total  energy  leads  to  a  reasonable  scaling 
for  Ec : 

*1 

Et  =  pe  +  —  p  (u2  +  v2) 

£t  =  pJ7i(p'e*  +  -|p*(  (u*)2  +  (V)2)) 

El  -  p’e*  +  -|p*(  (u*)2  +  (v*)2) 

=  E'tp’ul 

Note  that  in  the  freestream, 

El  =  ~  +  (  (Y-l)Y^-2)-1 

The  equations  of  motion  in  non-dimensional  form  turn  out  to  be 
the  same  as  the  dimensional  equations  when  the  asterisks  are 
dropped;  the  TVD  schemes  described  need  not  be  modified.  The 
Mach  number  only  enters  into  the  problem  through  the 
specification  of  Ec  in  the  far  field. 
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